function [T,Z]=rk4sgral(a,b,Za,M,fname)

% Método de Runge-Kutta para sistemas de ecuaciones 
% 
% Entrada
%         - a , b extremos del intervalo 
%         - Za=[x1(a), ..., xn(a)] condiciones iniciales
%         - M es el número de subintervalo en (a,b)
% Salida  - T el vector de valores equiespaciados
%		de la variable independiente
%         - Z=[x1(t) . . . xn(t)] donde xk(t) es la 
%	    aproximacion a la solución para la k_esima
%  		función incógnita.
% La función f(t,x) deberá ser ingresada en la línea de comando
% de Octave.

h=(b-a)/M;
Z=zeros(M+1,length(Za));
T=a:h:b;
Z(1,:)=Za;

for j=1:M   
	k1=h*feval(fname, T(j), Z(j,:));   
	k2=h*feval(fname, T(j)+h/2, Z(j,:)+k1/2);   
	k3=h*feval(fname, T(j)+h/2, Z(j,:)+k2/2);   
	k4=h*feval(fname, T(j)+h, Z(j,:)+k3);   

	Z(j+1,:)=Z(j,:)+(k1+2*k2+2*k3+k4)/6;

end
end
